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Preface 


A  technique  for  interpreting  Time  Domain  Electromagnetic  (TEM)  Data  was  proposed  in  ’’Detecting 
Unexploded  Ordnance  With  Time  Domain  Electromagnetic  Induction”  (Pasion,  1999).  An  approximate 
forward  modelling  for  the  TEM  response  of  compact  metallic  objects,  an  inversion  for  recovering  model 
parameters,  and  relationships  between  model  parameters  and  a  target’s  physical  parameters  were  estab¬ 
lished.  These  findings  were  combined  to  form  an  algorithm  for  locating  and  determining  the  approximate 
shape  of  a  buried  target. 


A  study  to  evaluate  the  effectiveness  of  the  technique  when  applied  to  Geonics  EM63  Time  Domain 
Electromagnetic  Sensor  data  was  undertaken  as  part  of  the  research  project  ’’Locating  and  Characterizing 
Unexploded  Ordnance  Using  Time  Domain  Electromagnetic  Induction”  (41262-EV),  funded  by  the  U.S. 
Army  Engineer  Research  and  Development  Centre  (ERDC).  A  series  of  TEM  measurements  was  carried 
out  between  April  5  and  April  21,  2000  at  the  ERDC  UXO  test  site  in  Vicksburg,  Mississippi.  The  data 
were  then  analysed  at  the  University  of  British  Columbia  Geophysical  Inversion  Facility  in  Vancouver, 
British  Columbia,  Canada. 


The  finHings  of  the  study  are  reported  here.  The  report  is  divided  into  two  parts.  Part  I  reproduces 
the  paper  ”A  Discrimination  Algorithm  For  UXO  Using  Time  Domain  Electromagnetics”  (Pasion  and 
Oldenburg,  2001)  which  describes  the  UXO  discrimination  algorithm  and  contains  an  application  of  the 
algorithm  to  one  of  the  data  sets  collected  at  the  ERDC  UXO  test  site.  Part  II  contains  descriptions 
of  the  different  Geonics  EM63  TEM  data  sets  acquired  during  the  study.  The  section  outlines  the  tests 
verifying  the  applicability  of  the  two-dipole  model  and  concludes  with  the  discrimination  algorithm  being 
applied  to  a  number  of  UXO  targets. 
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Part  I 


A  Discrimination  Algorithm  For  UXO  Using 
Time  Domain  Electromagnetics 


1  Abstract 


An  assumption  is  made  that  the  Time  Domain  Electromagnetic  (TEM)  response  of  a  buried  axisymmetric 
metallic  object  can  be  modelled  as  the  sum  of  two  dipoles  centered  at  the  midpoint  of  the  body.  The 
strength  of  the  dipoles  depends  upon  the  relative  orientation  between  the  object  and  the  source  field,  and 
also  upon  the  shape  and  physical  properties  of  the  body.  Upon  termination  of  the  source  field,  each  dipole 
is  assumed  to  decay  as  A;  (t  +  a)~^  The  parameters  k,  a,  /?  and  7  depend  upon  the  conductivity, 

permeability,  size  and  shape  of  the  object,  and  these  can  be  extracted  from  the  measurements  by  using 
a  nonlinear  parametric  inversion  algorithm.  Investigations  carried  out  using  an  analytic  solution  for  a 
sphere  and  laboratory  measurements  of  steel  and  aluminum  rectangular  prisms,  suggest  the  following 
two-step  methodology:  (1)  The  value  of  /?  is  first  used  as  a  diagnostic  to  assess  whether  the  metallic 
object  is  non-magnetic  or  magnetic,  (2)  the  ratios  of  fci/fca  and  P1/P2  are  then  diagnostic  indicators  as 
to  whether  the  geometry  is  plate-like  (uninteresting)  or  rod-like  (a  high  candidate  for  being  a  UXO). 
Results  from  the  application  of  this  algorithm  to  a  TEM  field  data  set  acquired  at  the  United  States 
Army  Corps  of  Engineers  Research  and  Development  Centre  (ERDC)  UXO  Test  Site  have  successfully 
identified  a  UXO  to  be  magnetic  and  rod-like. 
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Chapter  1 


Abstract 
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Introduction 


An  explosive  ordnance  is  a  munition  that  is  either  launched  or  fired  with  the  intent  of  detonation  at  a 
specified  target.  An  unexploded  ordnance  (UXO)  is  an  explosive  ordnance  that,  due  to  some  malfunction, 
remains  undetonated.  As  a  result,  the  ordnance  can  be  found  at  the  ground  surface,  partially  buried, 
or  buried  at  a  depth  of  up  to  8  m  beneath  the  surface.  The  remediation  of  UXO-contaminated  land 
has  been  made  a  high  priority  by  the  United  States  Department  of  Defense  in  order  to  either  maintain 
safe  usage  for  continuing  military  operations  or  to  permit  land  transfer  to  the  private  sector.  Practical 
and  cost-effective  strategies  for  remediation  require  both  detection  of  possible  targets  and  the  ability  to 
discriminate  between  UXO  and  contaminating  scrap  metal. 


The  detection  of  buried  metallic  objects  can  be  accomplished  with  a  variety  of  geophysical  sensing  tech¬ 
niques.  Time  domain  electromagnetic  induction  (TEM)  surveys  have  been  successful  in  detecting  both 
ferrous  and  non-ferrous  metallic  objects  near  the  soil  surface,  and  are  a  mainstay  amongst  technologies 
currently  utilized  in  UXO  clearance  projects.  In  the  TEM  method  a  time  varying  magnetic  field  is  used 
to  illuminate  a  conducting  target.  This  primary  field  induces  surface  currents  on  the  target  which  then 
generate  a  secondary  magnetic  field  that  can  be  sensed  above  ground.  With  time,  the  surface  currents 
diffuse  inwards,  and  the  observed  secondary  field  consequently  decays.  The  rate  of  decay,  and  the  spatial 
behavior  of  the  secondary  field,  is  determined  by  the  target’s  conductivity,  magnetic  permeability,  shape, 
and  size. 


Identification  of  a  UXO  from  electromagnetic  sensor  data  remains  a  major  hurdle  in  reducing  the  high 
costs  of  remediation  projects.  It  has  been  reported  that  approximately  70%  of  remediation  costs  are  cur¬ 
rently  being  used  to  excavate  non-ordnance  items  (Butler  et  al,  1998).  The  development  of  discrimination 
algorithms  can  be  roughly  categorized  as  either  model-based  or  data-based.  Data-based  algorithms  are 
pattern  recognition  procedures  that  compare  a  library  of  catalogued  responses  from  various  UXO  items 
to  measured  responses  (for  example  Damarla  and  Kessler,  2000).  Model-based  algorithms  use  either 
an  exact  or  approximate  forward  modeling  algorithm  to  determine  a  set  of  model  parameters  needed  to 
replicate  the  measured  responses,  and  subsequently  relating  the  model  parameters  to  physical  parameters 
(Khadr  et  al,  1998).  One  such  model-based  technique  that  has  been  the  focus  of  much  recent  research  is 
the  determination  of  the  time  constants  of  the  TEM  response,  or  equivalently  the  poles  of  the  frequency 
domain  signal,  to  identify  the  buried  target  (Snyder  et  al,  1999;  Baum,  1997;  Collins  et  al,  1999).  A 
method  that  represents  a  hybrid  of  the  model-based  and  data-based  algorithm  is  under  development  at 
Blackhawk  Geometries  (Grimm,  2000).  In  that  approach,  a  spheroid  modeller,  working  jointly  with  a 
model-based  inversion  algorithm,  generates  a  library  of  model  parameters  which  can  then  be  operated 
upon  by  a  neural  network  classifier  for  comparison  with  parameters  derived  from  the  raw  sensor  signal. 

In  this  paper  we  present  a  model-based  TEM  data  interpretation  algorithm  which  estimates  the  basic 
shape  (rod-like  or  plate-like)  and  magnetic  character  (ferrous  or  non-ferrous)  of  a  buried  metallic  object. 
We  first  present  an  approximate  forward  model  that  represents  the  time  domain  response  of  a  metallic 
object  as  a  pair  of  perpendicular  dipoles  located  at  the  center  of  the  buried  target.  This  form  of  model 
was  suggested  to  us  in  a  personal  communication  from  J.D.  McNeill.  The  strengths  of  these  dipoles 
decay  with  time,  and  the  parameters  that  govern  the  time  decay  behavior  are  related  to  the  conductivity. 
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permeability,  shape,  and  size  of  the  buried  target.  The  parameterization  is  simple,  and  thus  convenient  to 
use  in  data  fitting  procedures.  We  next  describe  an  inversion  scheme  to  recover  the  model  parameters  from 
TEM  data.  Since  these  parameters  encapsulate  information  about  the  physical  attributes  of  the  target, 
we  can  attempt  to  use  them  to  determine  if  the  target  is  ferrous  and  if  the  geometry  is  rod-like  (most 
likely  a  UXO)  or  plate-like  (most  likely  a  non-ordnance  item).  Empirical  relationships  are  developed  that 
link  the  model  parameters  to  the  physical  parameters  of  the  target  and  these  relationships  form  the  basis 
of  our  algorithm.  We  conclude  with  the  application  of  the  algorithm  to  a  synthetic  data  set  contaminated 
with  noise,  and  field  data  sets  taken  over  a  buried  UXO  and  buried  metallic  scrap. 
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3  Development  of  Approximate  Forward  Modelling 


In  order  to  invert  measured  TEM  data  for  the  physical  parameters  of  the  target,  it  is  necessary  to  have 
a  forward  model  to  describe  the  TEM  response  for  a  buried  metallic  object.  We  can  restrict  our  search 
for  response  solutions  to  axi-symmetric  metallic  targets,  since  this  geometric  subset  adequately  describes 
all  UXO  and  the  majority  of  buried  metallic  scrap  encountered  in  a  remediation  survey.  Unfortunately 
analytic  expressions  for  the  time  domain  response  are  restricted  to  a  metallic  sphere,  and  even  an  ex¬ 
pression  for  a  permeable  and  conducting  non-spherical  axi-symmetric  body  is  not  available.  Numerical 
solutions  of  Maxwell’s  equations,  under  continual  development,  are  promising  (e.g.  Haber,  2000;  Garin, 
2000;  Hiptmair,  1998),  however,  the  computational  time  requirements  for  obtaining  a  solution  still  make 
them  impractical  for  use  as  part  of  a  rigorous  inversion  procedure.  Our  approach,  therefore,  is  to  use  an 
approximate  forward  modelling  that  can  adequately  reproduce  the  measured  electromagnetic  anomaly 
in  a  minimal  amount  of  time.  The  validity  of  this  reduced  modelling  still  needs  further  testing  but  the 
empirical  tests  carried  out  here  suggest  that  it  can  be  useful  in  practice. 

The  development  of  the  approximate  forward  modelling  is  presented  in  four  steps.  We  begin  with  the 
response  of  a  sphere,  so  that  the  magnetic  polarization  dyadic  M  is  introduced.  This  dyadic  is  then  altered 
so  that  it  is  applicable  to  an  axi-symmetric  body.  This  generates  the  ”  two-dipole”  model  mathematically. 
Next  we  introduce  a  parameterization  for  the  time  decays  of  each  of  the  two  dipoles  and  finally,  we 
combine  everything  to  generate  our  approximate  forward  modelling. 


Response  for  a  Spherical  Body 

Consider  a  permeable  and  conducting  sphere  of  radius  a  illuminated  by  a  uniform  primary  field  (Fig. 
1(a))  At  a  time  t  =  0  the  primary  field  is  terminated,  and  eddy  currents  are  induced  in  the  sphere;  they 
subsequently  decay  due  to  the  finite  conductivity  of  the  sphere.  The  secondary  field  B®  generated  by  the 
decaying  currents  is  dipolar: 

B*«)  =  J^n.(<).(3f»-i)  (1) 

where  m  (t)  is  the  dipole  moment  induced  at  the  center  of  the  sphere  at  time  r  is  the  distance  between 
the  observation  point  and  the  sphere  center,  f  is  the  unit  vector  pointing  from  the  sphere  center  to  the 
observation  point  P,  and  I  is  the  identity  dyadic.  The  dipole  moment  is 

in(f)  =  —  B*’L®(f)  (2) 

Mo 


where 


L®  (t)  -  6aVr  E  ^  2) 


(3) 


^  ^  and  fir  =  is  the  relative  permeability  (Kaufmann  and  Keller,  1985).  In  pneral 

the  magnetic  permeability  of  highly  permeable  materials  is  a  function  of  many  parameters,  including  the 
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strength  of  the  incident  magnetic  field,  temperature,  and  magnetic  history.  However,  calculated  TEM 
responses  assuming  a  constant  permeability  of  fir  =  150  for  steel  and  fir  =  1  for  aluminum  compared 
well  with  laboratory  TEM  measurements  of  steel  and  aluminum  targets  (Pasion,  1999).  Therefore  we 
feel  that  Eqs.  (1)  to  (3)  are  suitable  for  the  analysis  that  follows  here.  The  values  g*  are  roots  to  the 
transcendental  equation 


tsm  Qg 


{fir  1)  Qs 

Qs  +  (Mr  -  1)  ' 


(4) 


Equations  (1)  to  (4)  reveal  that  the  secondary  B-field  of  a  sphere  in  a  uniform  primary  field  is  equivalent 
to  the  B-field  of  a  single  magnetic  dipole  located  at  the  center  of  the  sphere  and  oriented  parallel  to  the 
primary  field.  For  convenience  we  write  the  relationship  between  the  induced  dipole  and  the  primary 
field  as  m  =  •  B**,  where  M  is  the  magnetic  polarizability  dyadic.  For  a  sphere. 


2it 


27r 


M  =  — (t)I  = 

fj'o  Mo 


{t)  0 

0  L®  (t) 

0  0 


0 

0 


(5) 


Baum  (1999)  details  the  characteristics  of  the  magnetic  polarizability  dyadic,  and  notes  that  the  triple 
degeneracy  of  the  magnetic  polarizability  dyadic  reflects  the  symmetry  of  the  sphere. 


The  sphere  solution  possesses  several  characteristics  that  we  retain  in  the  formulation  of  our  approximate 
solution  for  an  axi-symmetric  target.  Firstly,  the  secondary  field  due  to  the  induced  currents  generated 
in  a  sphere,  illuminated  by  a  uniform,  step-off  primary  field,  is  dipolar  at  all  points  outside  the  sphere. 
We  will  also  represent  the  secondary  field  for  more  general  shapes  as  a  dipolar  field  (Eq.  (1)).  A  dipolar 
field  approximation  is  reasonable  for  any  observation  point  far  enough  away  from  any  localized  current 
distribution  (Jackson,  1975),  and  it  has  been  reported  that  for  observation  points  greater  than  1  to  2 
times  the  target  length,  a  dipolar  field  assumption  is  adequate  (Casey  and  Baertlein,  1999  or  Grimm  et 
ah,  1997).  Indeed,  higher  order  multipoles  induced  in  a  target  will  decay  at  early  times  (Grimm  et  ah, 
1997). 

Secondly,  the  induced  dipole  moment  in  the  center  of  a  sphere  is  given  by  the  dyadic  product  M  •  B^. 
This  form  indicates  that  the  induced  dipole  is  proportional  to  the  projection  of  the  primary  field  along 
the  direction  of  the  induced  dipole.  The  components  of  M  scale  the  strengths  of  the  dipoles.  The 
magnetic  polarizability  dyadic,  in  the  case  of  the  sphere,  contains  the  function  L®  {t)  that  contains  all 
the  information  about  the  time  decay  of  the  sphere  and  it  depends  upon  the  material  properties,  shape, 
and  size  of  the  target.  Our  hypothesis  is  that  more  general  metallic  shapes  can  also  be  approximately 
modelled  with  an  induced  dipole  equal  to  the  dyadic  product  M  •  B^.  However,  choosing  the  right 
functional  form  of  M  will  be  crucial. 


Approximating  the  Magnetic  Polarizability  Dyadic  for  an  Axi-Symmetric 
Body 

Analytic  expressions  for  M  for  the  time  domain  response  of  a  permeable  and  conducting  non-spherical  axi- 
symmetric  body  are  not  available.  Therefore  we  base  our  form  of  M  on  the  magnetostatic  j)olarizability 
for  a  spheroid.  Recall  that  for  the  time  domain  response  of  a  sphere  the  structure  of  M  is  identical 
to  the  structure  of  the  magnetostatic  polarizability  dyadic  of  a  sphere.  The  analytic  solution  for  the 
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magnetostatic  response  of  a  magnetic  prolate  spheroid  is  equivalent  to  the  field  of  a  magnetic  dipole 
induced  at  the  spheroid  center  (Das  et  al,  1990): 


^spheroid  =  irn  +  m2 


=  fci  [{z'  •  BP)  z']  +  k2  [(y'  •  BP)  y'  +  (x'  •  BP)  x'] 


>2  0 

0  k2 
0  0 


0 

0 

fci 


•BP 


(6) 


where  ifci  and  k2  are  the  polarizability  constants,  which  are  functions  of  the  conductivity,  permeability, 
shape,  and  size.  Equation  (6)  reveals  that  the  total  induced  dipole  can  be  written  as  the  sum  of  two 
orthogonal  dipoles  mi  and  m2.  The  first  dipole  moment  mi  is  parallel  to  the  major  axis  (z'  in  fig.  1(b))  of 
the  spheroid,  and  its  strength  is  proportional  to  the  product  of  the  primary  field  along  that  direction  and 
the  polarizability  fci .  The  second  dipole  moment  is  perpendicular  to  the  major  axis,  and  its  magnitude 
is  proportional  to  the  component  of  the  primary  field  along  that  direction  and  the  polarizability  k2.  A 
consequence  of  Jki  and  k2  being  functions  of  the  spheroid’s  shape  and  size  is  that  the  orientation  of  the 
effective  dipole  will  not  be  solely  determined  by  the  direction  of  the  primary  field,  as  is  the  case  for  a 
sphere.  In  addition,  the  orientation  of  will  be  influenced  by  the  aspect  ratio  of  the  spheroid. 


The  polarization  dyadic  in  Eq.  (6)  suggests  a  magnetic  polarization  dyadic  for  the  TEM  problem  of  the 
form 


M  = 


X2  (t) 
0 
0 


0 

L2  (t) 

0 


0 

0 

Li  it) 


(7) 


where  we  have  simply  replaced  fcj  and  ^2  in  Eq.  (6)  with  the  dipole  decay  functions  L\  (t)  and  L2  (t). 
The  resultant  induced  dipole  moment  for  this  definition  of  the  magnetic  polarization  dyadic  is  then 


m  (t)  =  mi  (t)  +  m2  it)  /g, 

=  Li  it)  [(z'  •  BP)  z']  +  L2  it)  [(y' .  BP)  y'  +  (x'  •  BP)  x'] 

Therefore,  our  approximate  forward  model  represents  the  TEM  response  of  two  orthogonal  dipoles.  The 
first  dipole  is  parallel  to  the  symmetry  axis  of  the  target,  and  the  second  dipole  is  perpendicular  to 
the  symmetry  axis.  These  dipoles  decay  independently  according  to  the  decay  laws  Li  (t)  and  L2  it), 
respectively. 


By  choosing  the  appropriate  parameters,  this  ’two-dipole’  model  produces  TEM  responses  that  are  con¬ 
sistent  with  those  observed  field  measurements  of  UXO.  It  has  been  noted  that  the  shape  anomaly  of  the 
measured  response  for  UXO  changes  with  time  (Grimm  et  al.,  1997).  The  physical  phenomenon  that  gave 
rise  to  the  temporal  changes  in  shape  anomaly  was  explained  in  terms  of  the  nature  of  the  induced  eddy 
currents.  Eddy  currents  that  circulate  end-to-end  in  the  UXO  dominate  at  early  time  but  decay  away 
quickly,  while  eddy  currents  that  circulate  about  the  long  axis  extend  later  into  time.  This  observed  field 
behavior  can  be  duplicated  by  letting  the  two  orthogonal  dipoles  mi  (t)  and  m2  (t)  decay  independently 
of  each  other.  The  dipole  mi  (t)  is  parallel  to  the  long  axis  and  it  simulates  the  magnetic  fields  that 
arise  from  currents  circulating  about  the  axis.  The  dipole  m2  (t)  is  perpendicular  to  the  long  axis  and 
it  simulates  the  magnetic  fields  that  arise  from  currents  circulating  end-to-end.  By  assigning  a  different 
decay  characteristic  (governed  by  its  decay  parameters)  to  each  dipole,  the  relative  contribution  by  each 
dipole  to  the  secondary  field  can  vary  with  time. 


Approximating  the  Magnetic  Polarizability  Dyadic  for  an  Axi-Symmetric  Body 
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Time  Decay  Functions  Li{t)  and  L2{t) 

The  time  decay  for  a  sphere  is  determined  by  the  sum  of  exponentials.  This  result  generalizes  to  the 
case  of  a  conductive  body  of  arbitrary  size  and  shape  in  an  insulating  medium  illuminated  by  a  step-off 
primary  field  (Kaufman,  1994).  Thus  the  form  for  L{t)  should,  at  least,  be  able  to  duplicate  the  time 
decay  features  observed  for  the  sphere.  Plots  of  the  B-field  and  dB/dt  response  for  both  a  magnetically 
permeable  (e.g.  steel)  and  non-permeable  (e.g.  aluminum)  sphere  are  shown  in  Fig.  2. 

An  appropriate  form  of  the  decay  law  for  the  B-field  is 

L{t)  =  k  {t  +  a)~^  e~ .  (9) 

The  parameter  k  controls  the  magnitude  of  the  modelled  response.  The  three  parameters  a,  P,  and  7, 
control  the  duration  and  characteristics  of  the  three  different  stages  of  the  time  decay  curve.  The  duration 
of  the  relatively  flat  early  time  stage  is  proportional  to  the  parameter  a.  The  linear  decrease  of  response 
observed  during  the  intermediate  time  stage  is  determined  by  The  exponential  decay  characterizing 
the  late  time  stage  is  controlled  by  the  parameter  7.  Figure  2(a)  demonstrates  the  ability  to  reproduce 
the  secondary  B-field.  This  form  of  the  decay  law,  with  the  a  parameter  absent,  was  suggested  to  us  in 
a  personal  communication  from  J.D.  McNeill. 

The  time  derivative  dB/dt,  which  is  measured  directly  with  most  TEM  receivers  can  also  be  modelled 
with  Eq.  (9).  Figure  2(a)  includes  plots  of  the  dB/dt  curves  for  a  steel  and  an  aluminum  sphere.  The 
early  time  behaviour  for  the  non-permeable  sphere  follows  a  decay  and  so  these  curves  are  different 
from  those  of  B  in  Fig.  2(b).  Nevertheless  the  curves  are  still  represented  by  early  time  turn-overs,  and 
linpar  and  exponential  decays  that  can  be  accommodated  by  Eq.  (9).  The  suitability  is  demonstrated  by 
the  fit  between  the  laboratory  measured  response  and  a  predicted  response  obtained  by  evaluating  Eq. 
(9). 

In  the  following  section  we  generically  denote  the  TEM  response  as  ^  (r,t)  where  ^  can  be  the  magnetic 
field  or  its  time  derivative.  The  time  dependent  decay  of  ^  is  given  by  Eq.  (9). 


The  Approximate  Forward  Model 


With  the  above  background,  we  can  write  an  approximate  expression  for  the  secondary  field  response 
of  an  axi-symmetric  target.  First,  let  us  switch  from  the  body-fixed  (primed)  coordinate  system  to  a 
space-fixed  coordinate  system,  which  is  more  amenable  to  the  definitions  of  target  and  sensor  location  of 
a  typical  field  survey  (Fig.  3).  A  vector  v'  in  the  body-fixed  system  co-ordinate  system  is  related  to  a 
vector  V  in  the  space-fixed  co-ordinate  system  via  the  Euler  rotation  tensor  A  {ip,  9, 4>)  by  (Arfken,  1985) 

v'  =  A  V.  (10) 


Due  to  the  axial  symmetry  of  the  problem  ip  =  0,  and  the  Euler  rotation  tensor  can  be  written 


A  = 


cos  9  cos  <p 
—  sin^ 
sin  9  cos  <p 


cos  9  sin  <p 
cos<p 
sin  9  sin  <p 


—  sin0 
0 

cos^ 


(11) 


where  6  is  the  angle  between  the  symmetry  axis  of  the  target  (z'  in  Fig.  1(b))  and  the  vertical  axis  in 
the  space-fixed  coordinate  system  (z  in  Fig.  3),  and  (j)  is  the  angle  between  the  projection  of  z'  onto  the 
horizontal  plane  and  x. 
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Time  Decay  Functions  Li{t)  and  1^2 (t) 


The  approximate  forward  modelling  is  written  by  substituting  the  definition  of  the  induced  dipole  of 
Eq.  (8)  into  the  expression  for  a  dipole  field  (1),  and  carrying  out  the  dyadic  product.  Let  us  consider  a 
target  whose  center  is  located  at  R  in  the  space-fixed  co-ordinate  system.  The  secondary  response  4  (r,  t) 
measured  at  a  receiver/transmitter  location  r  and  at  a  time  t  after  the  termination  of  the  primary  field, 
is  then  the  sum  of  the  responses  of  the  two  orthogonal  dipoles: 


where 


and 


e(r,t)=^i  (r,t)+^2  (r,t) 

(12) 

-  4;^  |r  -  Kf  |r  -  Rl"  ) 

(13) 

mi  (t)  =Li  (t)  (z'-BP)z' 

(14) 

m2  (t)  =L2  it)  [(x'  •  RP)  x'  +  (y'  •  RP)  y'] 

(15) 

are  the  dipole  parallel  and  perpendicular  to  the  axis  of  symmetry.  The  unit  vectors  are  given  by  Eq.  (11). 


In  summary,  the  approximate  response  of  a  buried  metallic  object  given  by  Eq.  (12)  can  be  generated 
from  13  parameters  that  describe  the  object.  These  model  parameters  are  elements  of  the  model  vector 

m  =  [X,  Y,  Z,  (p,  6,  ki,  ai,  /?i,  Ji,k2,  aa,  P2,  T2]  ■  (1®) 

X  and  Y  denote  the  surface  projection  of  the  centroid  of  the  body,  and  Z  is  the  depth  of  the  object 
below  the  surface.  The  orientation  of  the  target  is  described  by  the  two  angles  0  and  <t>.  The  remaining 
parameters  describe  the  decay  characteristics  of  the  two  dipoles:  ki,  ai,  and  71  describe  the  dipole 
parallel  to  the  axis  of  symmetry  (mi),  and  fc2,  0:21  P2,  and  72  describe  the  dipole  perpendicular  to  the 
axis  of  symmetry  (m2).  Thus  the  inversion  for  the  model  m  will  immediately  give  estimates  of  target 
location  and  orientation.  Information  on  the  shape,  size,  and  material  parameters  of  the  target  may  later 
be  inferred  from  the  remaining  parameters. 


The  Approximate  Forward  Model 
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4  Non-Linear  Parameter  Estimation  Procedure 


In  this  paper  it  is  first  assumed  that  the  response  meaisured  in  a  survey  is  due  to  a  single  body,  and 
second,  that  the  response  of  this  single  body  can  be  accurately  modelled  with  Eq.  (12).  With  these 
hypotheses,  an  inversion  procedure  is  developed  that  utilizes  the  approximate  forward  model. 

The  forward  model  can  be  expressed  as 

dj  =  Fj[ml  i  =  l,2,3,...iV  (17) 

This  equation  expresses  the  mapping  of  the  model  vector  m  to  a  datum  dj  by  a  functional  Fj.  The 
forward  mapping  Fj  is  defined  by  Eq.  (12)  and  is  a  nonlinear  functional  of  the  13  model  parameters  given 
in  Eq.  (16).  In  the  inverse  problem,  these  parameters  are  retrieved  from  a  vector  of  observed  data 
by  minimizing  a  least-squares  objective  function.  Before  proceeding  to  the  details  of  the  inversion  there 
are  two  important  practical  aspects  to  be  introduced:  (1)  Selected  parameters  must  remain  positive  and 
(2)  Parameters  should  be  scaled  to  enhance  stability  in  the  iterative  process. 

In  the  approximate  forward  model  the  time  decay  parameters  oii,  and  '7^  {i  1>2)  are  defined 
as  positive.  In  the  inverse  problem  the  positivity  of  these  parameters  can  be  maintained  by  solving  the 
associated  square- variable  unconstrained  problem  (Gill  et  al.,  1981).  Following  this  formulation,  each  time 
decay  parameter  mi  is  replaced  by  the  squared  variable  Wi  such  that  A  second  transformation 

is  a  linear  scaling  that  ensures  that  each  component  of  m  is  of  order  unity.  This  is  done  by  dividing  each 
parameter  mt  by  its  typical  value.  In  the  following  presentation  of  the  inversion  algorithm,  the  letter  m 
continues  to  denote  the  parameter  vector  to  be  inverted  for,  but  now  includes  the  positivity  and  scaling 
transformations . 


Defining  the  Objective  Function 

If  there  axe  L  time  channels  and  K  locations  where  TEM  data  are  collected,  then  there  will  be  iV  =  KL 
data  points  contained  in  the  data  vector  d®'’*.  Because  data  will  be  collected  on  several  lines,  with  a 
number  of  stations  per  line,  there  will  generally  be  far  more  data  than  model  parameters  (N  »  13). 
Therefore  the  inversion  for  m  involves  solving  an  overdetermined  system  of  non-linear  equations,  with 
the  goal  of  finding  the  model  that  produces  the  data  that  best  fit  the  observed  data.  This  is  a  non-linear 
least  squares  problem  and  is  solved  by  minimizing 

#(m)  =  i||W<,(F[m]-d°^*)|F  (18) 

where  E[m]  are  the  forward  modelled  data,  d®'’*  are  the  observed  data,  and  #  is  the  least  squares  objective 
function  that  measures  how  closely  our  predicted  data  match  the  observed  data.  Wd  is  the  data  weighting 
matrix.  If  the  data  are  contaminated  with  unbiased  Gaussian  random  noise,  then  Wj  is  ideally  a  diagonal 
matrix  whose  elements  are  the  reciprocals  of  the  standard  deviation  of  each  datum.  The  noise  arises 
from  many  sources,  including  sensor  location  errors,  instrument  noise,  and  inaccuracy  of  the  forward 
modelling.  It  is  likely  that  the  Gaussian  independent  assumption  is  not  valid,  but  it  is  essential  to 
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estimate  a  quantity  that  reflects  the  uncertainty  in  each  datum.  We  assume  that  the  errors  can  be 
characterized  by  a  percentage  of  the  datum  value  plus  a  threshold,  that  is 


(Wd),,  = 


1 

pdi  +  e 


(19) 


where  p  is  typically  a  percentage  and  e  is  a  constant  that  characterizes  ambient  noise.  The  positive  e 
ensures  that  small  data  points  would  have  reasonable  errors  assigned  to  them,  and  thus  prevents  them 
from  having  undue  influence  on  the  solution. 


Minimizing  the  Objective  Function 

For  ease  of  notation  the  least  squares  problem  is  rewritten  as 

1  1  ^ 

minimize  $  (m)  =  -r  (m)^  r  (m)  =  -  ^  ri  (m)^  (20) 

2  ^  i=i 

where  r  is  the  residual  function 

r  (m)  =  Wd  (F[m]  -  d'’**) 

and  ri  (m)  is  the  component  of  r  (m).  We  adopt  a  modified  Newton’s  method  to  minimize  the  objective 
function.  The  approach  taken  here  is  to  first  make  an  initial  guess  of  the  model  parameters  m^.  Techniques 
for  making  this  initial  guess  are  outlined  in  Pasion  (1999).  The  starting  model  is  iteratively  improved 
to  find  the  minimum  of  Eq.  (20).  At  each  iteration  a  Newton’s  search  direction  is  chosen  that  minimizes 
the  local  quadratic  model  about  the  current  iterate  mj^.  The  Newton  step  Sm  for  the  non-linear  least 
squares  problem  is  then  given  by 

H(mfc)^m  = -J(mA)^r(m*)  (21) 

where  the  Jacobian  matrix  J  and  the  Hessian  matrix  H  are  defined  as 
Jy  (m)  =  and  H  =  J  (m)’’  J  (m)  +  S  (m) , 

C/XT\j 

where  S  (m)  is 


N 

s  (m)  =  ^  n  (m)  V^r,  (m) 

i=l 

A  new  model  =  m*  +  X6m  is  then  defined,  where  the  positive  scalar  A  is  chosen  such  that 

$  +  XSxn)  <  ^  (rrijt).  The  sequence  of  iterations  is  terminated  once  the  relative  gradient  measure 

is  less  than  a  tolerance  level  ,  or  once  there  is  insignificant  change  in  the  models  obtained  at  successive 
iterations  (Dennis  and  Schnabel,  1983). 


Error  Bounds  of  the  Parameter  Estimates 

Once  the  model  parameters  m,  which  minimize  the  objective  function  $  (m)  have  been  obtained,  we  can 
examine  the  reliability  and  precision  of  the  estimated  parameters  via  the  model  covariance  matrix.  Let 
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m+  be  the  best  estimate  of  the  model  in  the  absence  of  noise,  and  (5m*  —  m+  —  m*.  The  model  covariance 
matrix  Vm  is  defined  as  the  expectation  value  of  (5m*  (5m*^(Bard,  1974) 

y,n  =  E  ((5m*  (5m*’’)  «  E  (h*"'  J*’’  (5d'’'*  (5d‘’'’*’’j*H*-') 


The  Hessian  and  Jacobian  in  the  above  expression  are  evaluated  at  m  —  m*,  and  are  therefore  constants. 
As  a  result  they  can  be  taken  outside  of  the  expectation  value  expression: 


=  H*-iJ*’' Vrf  J*H*-' 


(22) 


where  is  the  covariance  matrix  of  the  data.  In  the  case  when  the  observations  have  uncorrelated  errors, 
the  data  covariance  matrix  reduces  to  a  diagonal  matrix,  and  an  estimate  of  the  standard  deviation  of 
the  f*'*  model  parameter  rrii  is  then 


=  (V„)ii  = 


J*^J,H* 


■) 


ii 


(23) 


Model  variance  estimates  applied  to  non-linear  problems  are  not  as  reliable  as  when  implemented  in  linear 
least  squares  problems,  and  they  should  only  be  used  as  a  very  rough  estimate  (Bard,  1974;  Dennis  and 
Schnabel,  1983).  Nevertheless,  Eq.  (23)  at  least  provides  a  minimum  estimate  to  the  uncertainties  of  the 
parameters. 
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5  Relating  Model  Parameters  to  Material  and 
Geometric  Properties 


The  above  inversion  generates  the  parameters  that  characterize  a  target’s  TEM  anomaly.  The  next  step 
is  to  interpret  these  parameters.  Recall  that  UXO  are  typically  rod-like  rather  than  plate-like,  and  are 
magnetically  permeable.  In  order  to  extract  these  potentially  UXO  identifying  features  from  the  recovered 
model  m*,  we  use  the  inversion  procedure  to  fit  a  series  of  decay  curves  from  a  range  of  axi-symmetric 
targets  of  different  shape,  geometry,  and  material  properties.  We  then  generate  empirical  relationships 
between  the  parameters  and  target  characteristics.  The  data  curves  used  for  this  analysis  were  either 
TEM  measurements  made  in  the  Geonics  Ltd.  laboratory,  or  they  were  synthetically  generated  decay 
curves  for  a  sphere  using  Eq.  (12). 


Lab  Setup  and  Measurements 

A  series  of  TEM  measurements  of  metallic  targets  was  made  by  Geonics  Ltd.  A  40  m  x  40  m  square 
transmitter  loop  was  used  to  provide  a  relatively  uniform  field  at  the  center  of  the  loop.  Aim  diameter 
receiver  coil  was  placed  coaxial  and  coplanar  to  the  transmitter  loop,  and  each  target  was  located  at  the 
center  of  the  receiver  loop.  The  Geonics  PROTEM  47  time  domain  equipment  was  used  for  producing 
the  transmitting  field  and  for  recording  the  time  domain  measurement  due  to  a  step-off  current.  Mea¬ 
surements  of  the  time  decay  response  of  these  targets  were  recorded  as  plots  of  log(9B/5t)  vs.  log(t). 
Since  values  were  not  recorded  by  a  data  logger,  the  plots  were  subsequently  digitized.  Plots  of  the 
steel  target  responses  were  digitized  by  hand  by  J.D.  McNeill  at  Geonics  Ltd.,  and  the  aluminum  target 
responses  were  digitized  at  UBC  after  scanning  the  plots  into  a  computer.  Analyses  were  performed 
on  both  the  impulse  (9B/5t)  response  measured  as  an  induced  voltage  in  the  receiver,  and  also  on  the 
B-field  response.  The  B-field  response  was  obtained  by  integrating  the  induced  voltage. 

Two  sets  of  targets  were  measured.  The  first  set  of  measurements  involved  recording  the  TEM  response 
for  a  series  of  steel  and  aluminum  rectangular  prisms  of  different  dimensions.  Each  prism  had  at  least 
one  dimension  of  20.3  cm,  and  the  targets  ranged  from  a  thin  rod  (20.3  x  0.6  x  0.6  cm)  to  a  cube 
(20.3  X  20.3  X  20.3  cm)  to  a  thin  plate  (20.3  x  20.3  x  0.6  cm).  A  second  set  of  measurements  was  made 
on  24  sample  UXO.  These  targets  included  various  ordnance  items  used  by  NATO  since  World  War  II. 
The  ordnance  range  in  length  from  18  to  85  cm,  and  in  diameter  from  6.1  to  15.9  cm.  A  diagram  of  aU 
the  ordnance,  along  with  a  table  listing  the  dimensions  of  each  ordnance,  is  included  in  Pasion  (1999). 

The  axi-symmetric  targets  were  placed  in  two  orientations  at  the  center  of  the  receiver  loop  in  two 
orientations.  Each  target  was  measured  with  the  axis  of  symmetry  perpendicular  and  parallel  to  the 
primary  field.  Since  the  strength  of  each  induced  dipole  is  proportional  to  the  projection  of  the  primary 
field  onto  the  dipole  direction,  the  two  measurement  orientations  isolate  the  decay  behaviour  of  each  of 
the  two  dipoles.  For  example,  consider  a  plate.  When  the  primary  field  is  perpendicular  to  the  plane 
of  the  plate,  the  projection  of  the  primary  field  onto  dipole  2  is  zero,  thus  the  approximate  forward 
model  assumes  the  response  can  be  modelled  as  a  single  dipole  perpendicular  to  the  plate.  The  decay 
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parameters  of  dipole  1  (fci,  qi,  /3i,  and  71)  can  then  be  estimated  by  fitting  this  curve  to  the  decay  law 
(Eq.  (9)).  When  the  primary  field  is  parallel  to  the  plane  of  the  plate,  the  response  is  due  to  dipole  2 
and  parameters  k^,  0L2,  P2,  and  72  can  be  recovered. 


Relationships  between  the  target  characteristics  and  the  model  parameters  were  established  in  the  fol¬ 
lowing  manner.  A  scaled-down  version  of  the  non-linear  least  squares  techniques  outlined  in  the  previous 
section  was  used  to  obtain  the  decay  parameters  k,  a,  P,  and  7  for  each  of  the  target’s  two  dipoles.  Sec¬ 
ondly,  we  observed  how  recovered  values  of  model  parameters,  or  combinations  of  parameters,  changed 
with  the  dimensions  and  magnetic  properties  of  the  measured  prism.  The  patterns  in  the  behaviour  of 
the  parameters  then  led  to  the  shape  and  permeability  discrimination  diagnostics  that  are  proposed  in 
the  following  sections. 


A  Relationship  Between  ^  and  Magnetic  Permeability 

UXO  are  generally  made  of  steel,  which  is  a  ferrous  material.  Therefore,  the  magnetic  permeability  is 
likely  to  be  an  identifying  characteristic  of  UXO.  To  generate  a  link  between  magnetic  permeability  and 
model  parameters,  forward  modelled  responses  were  calculated  for  a  series  of  spheres  varying  in  size  and 
permeability.  Both  B-field  and  dB/dt  data  were  then  inverted  to  generate  decay  parameters,  and  in 
particular,  to  produce  estimates  of  the  parameter  p.  The  plots  of  /3  as  a  function  of  sphere  radii  and 
magnetic  permeability,  are  provided  in  Fig.  4. 

Figure  4  suggests  that  the  value  of  P  obtained  for  a  sphere  may  be  diagnostic  in  determining  whether 
the  sphere  is  permeable  or  non-permeable.  Figure  4(a)  exhibits  the  relationship  for  the  dB/ dt  responses. 
For  a  sphere  with  fir  =  150  (typical  for  steel),  we  see,  for  spheres  with  radii  between  5  to  15  cm,  that  P 
falls  between  1.11  and  1.35;  while  for  a  sphere  with  =  1,  has  a  value  of  approximately  0.5,  which 
corresponds  to  the  early  time  behavior  that  Kauhnan  (1994)  predicted  for  a  non-permeable  sphere. 
Therefore,  when  applying  our  inversion  to  the  time  derivative  of  the  field,  a  value  of  p^^!^  greater  than 
about  0.8  indicates  that  the  target  is  most  likely  permeable.  This  analysis  is  repeated  on  the  forward 
modelled  B-field  responses,  and  the  results  are  plotted  in  Fig.  4(b).  A  threshold  value  of  0^  =  0.3  could 
be  used  such  that  targets  with  a  B-field  response  characterized  hy  a,  P  >  0.3  indicates  a  permeable  target. 

The  use  of  /S  as  a  diagnostic  to  determine  permeability  can  be  extended  to  non-spherical  targets  by 
looking  at  the  recovered  P  values  for  the  aluminum  and  steel  prisms.  The  inversion  produces  two  values 
of  P,  one  for  each  of  the  excited  dipoles,  to  describe  a  buried  target.  We  suggest  taking  the  average  of  the 
two  recovered  P  values,  which  we  label  as  p.  When  analyzing  the  dBjdt  responses  of  the  axi-symmetric 
aluminum  targets,  =  0.52  with  a  standard  deviation  of  0.07.  For  the  steel  targets  =  1.11 

with  a  standard  deviation  of  0.08.  These  averages  fall  on  either  side  of  the  0.8  threshold  obtained  by 
fitting  sphere  dBjdt  responses. 

When  analyzing  the  B-field  responses  of  the  axi-symmetric  aluminum  targets,  P^  =  0.17  with  a  standard 
deviation  of  0.03.  For  the  steel  targets  P^  =  0.5  with  a  standard  deviation  of  0.2.  These  averages  fall  on 
either  side  of  the  0.3  threshold  obtained  by  fitting  sphere  B-field  responses,  and  so  again,  a  consistent 
criterion  can  be  used. 
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A  Relationship  Between  P  and  Magnetic  Permeability 


Relationships  Between  Model  Parameter  Ratios  and  Target  Shape 

Empirical  relationships  were  also  established  between  the  target  shape  and  the  ratios  fei/fc2  and  ^i//?2- 
For  space  reasons,  we  present  only  the  analysis  of  the  response  of  the  targets,  and  refer  the 

reader  to  Pasion  (1999)  for  the  analysis  of  the  B-field  data.  There  we  show  that  the  same  target  shape 
diagnostics,  developed  here  for  9B/9t  data,  also  apply  to  B-field  data. 


The  Ratio  kx/k^ 

The  recovered  k  values  for  targets  ranging  from  a  steel  plate  to  a  steel  rod  are  shown  in  Fig.  5(a),  and 
the  calculated  /k-ratios  are  shown  in  Fig.  5(b).  For  a  steel  plate,  the  fe-ratio  fci/A:2  <  1-  For  a  steel  bar 
the  jfc-ratio  fei/fe  >  1-  The  recovered  k  values  for  aluminum  targets  are  shown  in  Fig.  5(c).  The  opposite 
orientation  effect  was  observed  for  an  aluminum  rod,  that  is  fci/fc2  <  1  (Fig.  5(d)). 


The  Ratio  ^y/ ^2 

In  addition  to  the  relative  strength  of  the  dipoles  being  shape  dependent,  the  slope  of  the  time  decay 
response  (either  dB/dt  or  B-field)  during  the  intermediate  time  stage  is  dependent  upon  the  target  shape. 
This  effect  was  seen  in  steel  targets  only.  The  steepness  of  the  response  during  the  intermediate  time 
stage  is  reflected  in  the  parameter  yd.  The  recovered  yd  values  for  targets  ranging  from  a  steel  plate  to 
a  steel  rod  are  shown  in  Fig.  6(a),  and  yd  values  for  aluminum  targets  are  shown  in  Fig.  6(c).  A  dipole 
that  decays  at  a  greater  rate  will  have  a  larger  yd.  The  rate  of  decay  of  the  dBjdt  response  is  greater 
when  the  plane  of  a  steel  plate  is  perpendicular  to  the  primary  field  (dipole  1),  than  when  the  plane  of  a 
steel  plate  is  parallel  to  the  primary  field  (dipole  2).  Thus,  for  a  steel  plate  the  ^-ratio  ^i//?2  >  !•  In  th© 
case  of  a  rod,  the  dB/dt  response  decays  faster  (and  thus  yd  is  larger)  when  the  main  axis  of  the  rod  is 
perpendicular  to  the  primary  field  (dipole  2).  In  the  case  of  a  steel  rod  the  ^-ratio  ^1/^2  <  1  (Fig.  6(b)). 

For  aluminum  targets  the  response  shape  looks  essentially  the  same  for  each  of  the  targets.  The  dB/dt 
response  exhibits  a  power  law  decay  of  and  is  exponential  at  later  times.  The  decay  curves  for 

aluminum  targets  are  essentially  the  same  regardless  of  target  shape,  and  therefore  there  is  no  relationship 
between  the  yd-ratio  and  the  aspect  ratio  (Fig.  6(d)). 


The  Discrimination  Algorithm  Using  d^/dt  Data 


The  results  from  the  previous  section  suggest  the  following  algorithm  for  using  dB/ dt  data  to  help  identify 
possible  UXO  targets: 


1.  Perform  the  non-linear  inversion  outlined  in  the  previous  section  to  recover  model  parameters  for 
the  two-dipole  model. 

2.  Compute  yd  =  |(ydi  +  0.80  then  the  target  is  most  likely  permeable. 

3.  Compute  ratios  ydi/yd2  and  ki/k2.  There  are  two  options: 

•  yd  >  0.8  =>■  Ferrous  Target.  \i  ky/k2  >  1  and  ydi/yd2  <  1  then  a  permeable  rod-like  target  was 
measured.  If  ki/k2  <  1  and  ydi/^  >  1  then  a  permeable  plate-like  target  was  measured. 


Relationships  Between  Model  Parameter  Ratios  and  Target  Shape 
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,0  <  0.8  Non-Ferrous  Target:  If  fci/fca  >  1  then  non-permeable  plate-like  target  was  mea¬ 
sured.  If  ki/k2  <  1  then  the  target  is  rod-like.  P1/P2  does  not  give  supporting,  or  extra, 
information. 


The  above  algorithm  can  be  extended  to  the  analysis  of  B-field  data  simply  by  changing  the  0  threshold 
to  0.3. 

The  parameters  a  and  7,  which  respectively  are  representative  of  a  target’s  early  and  late  time  charac¬ 
teristics,  are  not  used  in  the  discrimination  algorithm.  The  a  parameter  is  very  sensitive  to  how  early  the 
measurement  of  the  TEM  sensor  extends  in  time  and  we  have  found  that  it  is  not  as  robustly  estimated 
as  fc,  and  7.  The  late  time  constant  for  conductors  is  related  to  a  target’s  shape,  size  and  conductivity 
(Kaufman,  1994  and  Nabighian  and  Macnae,  1991).  However,  as  noted  by  Nabighian  and  Macnae  (1991), 
using  the  time  constant  for  discrimination  is  difficult  without  prior  knowledge  of  the  shape  and  either  a 
length  scale  or  conductivity.  Our  tests  thus  far  have  been  unable  to  generate  a  consistent  relationship 
between  the  time  constant  recovered  from  lab  TEM  measurements  and  target  characteristics. 
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The  Discrimination  Algorithm  Using  dB/dt  Data 


6  Synthetic  Data  Set 


The  parameter  estimation  procedure  is  now  tested  on  a  synthetically  generated  field  data  set.  The 
object  of  interest  is  a  75  mm  anti-tank  mortar.  The  primary  decay  curves  for  the  axial  and  perpendicular 
orientations  were  obtained  from  measurements  made  by  Geonics  using  the  setup  described  in  the  previous 
section.  These  were  inverted  to  recover  the  decay  parameters  for  each  dipole.  For  this  simulation,  the 
target  is  assumed  to  be  buried  at  a  depth  of  57  cm  (Z=l  m),  and  located  at  (2m  N,2m  E)  on  the  survey 
grid.  The  mortar  is  oriented  such  that  =  30°  and  0  =  65°.  Equation  (12)  was  used  to  generate  the 
“observed”  data  for  this  example. 

The  survey  consists  of  a  2  m  x  2  m  grid,  containing  5  lines  of  data  separated  at  50  cm  line  spacing,  with 
stations  located  at  20  cm  intervals  along  each  line.  At  each  station  the  vertical  component  of  the  voltage 
is  generated  for  26  logarithmically  spaced  time  channels.  The  time  channels  range  from  0.01  ms  to  100 
ms.  In  order  to  make  this  example  closer  to  a  real  TEM  data  set,  5%  random  Gaussian  noise  was  added 
to  the  data  and,  since  a  real  TEM  instrument  will  have  a  finite  measurement  sensitivity  to  the  secondary 
field,  a  data  threshold  of  0.001  is  set.  There  are  1278  total  data  points  exceeding  the  minimum  threshold. 


The  inversion  is  carried  out  with  a  data  weighting  matrix  in  Eq.  (19)  with  p  =  0.05  and  e  =  0.001  mV. 
The  observed  data,  and  data  predicted  by  the  recovered  model,  are  compared  in  Figs.  7  and  8.  Figure  8 
shows  a  plan  view  comparison  for  three  of  the  26  time  channels.  The  difference  maps  exhibit  a  random 
distribution  over  the  data,  indicating  the  reluctance  of  the  inversion  to  fit  the  noisy  portion  of  the  data.  A 
comparison  of  the  true  model  iritrue,  recovered  model  m^ec  and  the  estimated  model  standard  deviations 
ct"*  are  found  in  Table  1.  The  discriinination  algorithm,  when  applied  to  the  recovered  decay  parameters, 
yields  the  following.  The  value  of  =  1.07(>  0.8)  indicates  that  the  target  is  likely  permeable.  The 
ratios  ki/k2  =  3.71(>  1)  and  ^i//32  =  0.71(<  1)  indicate,  for  a  magnetically  permeable  target,  that  the 
TEM  response  is  likely  from  a  rod-like  target. 


Chapter  6  Synthetic  Data  Set 
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7  Field  Data  Set 


We  now  apply  our  algorithm  to  a  TEM  field  data  set  acquired  at  the  ERDC  UXO  test  site  in  Vicksburg, 
Mississippi.  The  Geonics  EM63  instrument  used  for  the  survey  is  a  multi-time  channel  time  domain  unit 
consisting  of  a  1  m  x  1  m  square  transmitter  coil  and  a  coaxial  horizontal  circular  receiver  loop  mounted 
on  a  two-wheel  trailer.  Measured  voltages  are  averaged  over  26  geometrically  spaced  time  gates,  spanning 
the  range  0.18  ms  to  25.14  ms. 


A  105  mm  projectile  is  placed  in  the  ground  with  its  center  at  2.0  m  East,  1.83  m  North  and  at  a  depth 
of  0.44  m  from  the  surface.  The  projectile  was  placed  horizontal  (0  =  90°),  with  its  tip  pointing  to 
the  North  (</>  =  0°).  Once  the  target  was  placed  in  the  ground,  it  was  not  covered  in  soil.  The  survey 
consisted  of  a  2  m  x  2  m  grid  centered  on  the  target,  containing  5  lines  running  North-South  separated 
at  50  cm  line  spacing,  with  stations  located  at  5  cm  intervals  along  each  line.  A  measured  signal  of  less 
than  1  mV  is  assumed  to  be  indistinguishable  from  the  noise.  The  resulting  data  set  contains  1882  total 
data  points. 


The  inversion  is  carried  out  with  a  data  weighting  matrix  in  Eq.  (19)  with  p  =  0.05  and  e  —  1  mV.  The 
first  stage  of  the  time  decay  evident  in  Fig.  2  is  not  observed  in  the  time  window  recorded  by  the  EM63. 
Therefore,  we  invert  these  data  by  setting  o  to  be  a  small  constant  that  does  not  affect  the  predicted 
data  within  the  EM63  time  range.  So  only  parameters  fc,  P,  and  7  for  each  dipole  are  recovered.  The 
observed  and  predicted  data  are  compared  in  Figs.  9  and  10.  Figure  10  shows  a  plan  view  comparison  for 
five  of  the  26  time  channels.  At  early  times  the  anomaly  has  a  single  peak  located  approximately  above 
the  UXO  center.  This  peak  splits  into  two  distinct  peaks  at  late  time.  The  recovered  model  predicts  data 
that  reflects  this  behavior.  Figure  9  compares  the  decay  curve  measured  at  four  stations  on  the  survey. 

The  recovered  location  and  orientation  parameters  are  listed  in  Table  2(a).  The  recovered  easting  of  2.04 
m  differs  from  the  true  value  of  2.00  m  by  4  cm.  The  recovered  northing  of  1.77  m  differs  from  the  true 
value  of  1.83  cm  by  6  cm,  placing  the  inducing  dipole  closer  to  the  projectile  tail.  These  errors  are  of 
the  same  magnitude  as  can  be  expected  in  spotting  the  station  location  in  the  field  survey.  In  addition, 
the  buried  105  mm  projectile  has  a  copper  rotating  band  near  the  tail  of  the  projectile.  It  has  been 
suggested  that  the  presence  of  the  rotating  band  will  shift  the  location  of  the  induced  dipole  from  the 
target  center  towards  the  tail  (Miller,  2000).  The  recovered  burial  depth  of  0.47  m  is  3  cm  deeper  than 
the  expected  depths  of  0.44  m.  The  orientation  parameters  6  and  cj)  are  well  recovered.  The  recovered 
decay  parameters  are  listed  in  Table  2(b)  and  the  diagnostics  applied  to  these  parameters  are  listed  in 
Table  2(c).  The  value  of  ^  =  0.91(>  0.8)  indicates  that  the  target  is  likely  to  be  magnetically  permeable. 
The  ratios  ki/k^  =  2.63(>  1)  and  P1IP2  =  0.69(<  1)  indicate,  for  a  magnetically  permeable  target,  that 
the  TEM  response  is  likely  to  be  from  a  rod-like  target. 
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Chapter  7  Field  Data  Set 


8  Discussion  and  Conclusion 


Efficient  remediation  of  areas  containing  UXO  first  requires  that  purely  conductive  metal  targets  be 
distinguished  from  steel  targets  that  are  conductive  and  permeable.  The  second  stage  then  focuses 
upon  determining  if  the  steel  target  is  rod-like  (and  likely  to  be  a  UXO)  or  plate-like  (and  not  be  of 
interest).  To  attack  this  problem  we  propose  a  modified  parametric  model  from  which  TEM  responses 
can  be  estimated.  The  TEM  response  of  a  buried  axisymmetric  metallic  object  is  modelled  as  the  sum 
of  two  dipoles  located  at  the  midpoint  of  the  body.  Non-linear  inversion  methods  are  used  to  extract  the 
parameters  from  the  field  data  and  these  parameters  are  subsequently  used  in  a  discrimination  procedure 
which  has  two  parts.  First,  the  decision  about  whether  the  object  is  ferrous  might  be  made  by  examining 
the  size  of  the  recovered  /3  values.  Second,  if  the  object  is  considered  to  be  ferrous,  then  the  ratios  of  ki/k2 
and  /0i/(02  are  diagnostic  indicators  of  whether  the  geometry  is  plate-like  or  rod-like.  These  diagnostics 
were  developed  for  both  B-field  and  dJifdt  TEM  data. 

This  algorithm  was  applied  to  a  synthetic  data  set  as  well  as  to  a  field  data  set  collected  by  the  Geonics 
EM63  time  domain  electromagnetic  sensor  over  a  105  mm  projectile.  In  both  cases  the  diagnostics, 
applied  to  the  recovered  model  parameters,  correctly  predicted  that  the  TEM  anomaly  was  produced 
by  a  magnetically  permeable  and  rod-like  metallic  target.  Although  further  testing  will  be  required  to 
fully  evaluate  our  proposed  technique,  the  results  presented  here  are  promising  and  may  have  a  positive 
impact  on  the  interpretation  of  UXO  detection  data. 


Chapter  8  Discussion  and  Conclusion 
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Figure  1:  The  body-fixed  (primed)  coordinate  system  for  a  sphere  and  a  spheroid. 


Figure  2:  (a)  The  time  decay  behaviour  of  the  time  derivative  of  the  magnetic  field  d^/dt.  (b)  The  time  decay 
behaviour  of  the  magnetic  flux  density  B.  The  B-field  response  is  normalized  by  the  strength  of  the  primary 
field.  The  solid  lines  are  responses  evaluated  from  eq.  (9).  The  agreement  supports  the  validity  of  this  parametric 
representation  of  the  time  domain  responses. 
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Chapter  8  Discussion  and  Conclusion 


Figure  3:  The  field  (unprimed)  co-ordinate  system  for  a  buried  target.  The  unit  vectors  x,  y,  and  z  define  the 
field  co-ordinate  system,  and  x',  y',  and  z'  define  the  body-fixed  co-ordinate  system. 


Figure  4:  The  behaviour  of  parameter  ^  for  various  size  spheres  with  varying  permeability  p.  Panel  (a)  contains 
results  of  recovering  by  fitting  the  dB/dt  data.  Panel  (b)  contains  results  of  recovering  by  fitting  the  B-field. 
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Figure  5:  Relating  the  aspect  ratio  of  a  steel  target  with  the  ratio  ki/kz.  Plot  (a)  contains  the  recovered  k 
parameter  from  fitting  the  measured  dB/dt  response  of  steel  axi-symmetric  targets.  Plot  (b)  illustrates  the 
relationship  between  the  ki/ki  ratio  derived  from  dB/dt  data  and  the  shape  of  a  steel  target.  Plot  (c)  contains 
the  recovered  k  parameter  from  fitting  the  measured  dB/dt  response  of  aluminum  axi-symmetric  targets.  Plot 
(d)  illustrates  the  relationship  between  the  ki/k2  ratio  and  the  shape  of  an  aluminum  target. 
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Chapter  8  Discussion  and  Conclusion 


Figure  6:  Relating  the  aspect  ratio  of  a  steel  target  with  the  ratio  0x1  Plot  (a)  contains  the  recovered  P 
parameter  from  fitting  the  measured  response  of  steel  axi-symmetric  targets.  Plot  (b)  illustrates  the 

relationship  between  the  Pi /Pi  ratio  derived  from  dBldt  data  and  the  shape  of  a  steel  target.  Plot  (c)  contains 
the  recovered  p  parameter  from  fitting  the  measured  d^/dt  response  of  aluminum  axi-symmetric  targets.  Plot 
(d)  illustrates  the  relationship  between  the  pi  / p2  ratio  and  the  shape  of  an  aluminum  target. 
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(a)  Location  and  Orientation  (b)  Decay  Parameters 


rrii 

^true 

VTirec 

cr'" 

ki 

7.07 

12.02 

12.064 

0.18 

ax 

0.01 

0.0076 

0.00759 

0.0030 

pi 

1.00 

0.89 

0.890 

0.0070 

7i 

3.16 

17.65 

17.635 

0.16 

k2 

7.07 

3.30 

3.252 

0.071 

Ot2 

0.01 

0.0077 

0.0076 

0.0027 

P2 

1.00 

1.25 

1.252 

0.014 

72 

3.16 

11.54 

11.68 

0.33 

rrii 

mf.j'ue 

nirec 

Northing  (m) 

1.90 

2.00 

2.00060 

0.00254 

Easting  (m) 

2.15 

2.00 

2.00105 

0.00169 

Depth  from 
loop  (m) 

1.20 

1.00 

1.00051 

0.00764 

(j)  (degrees) 

45 

30 

30.07 

0.18 

6  (degrees) 

45 

65 

65.04 

0.11 

Table  1:  Recovered  parameters  from  the  inversion  of  the  synthetic  data  set,  m©  is  the  starting  model.  The  true 
model  mtrue  and  the  recovered  model  rrirec  are  close.  The  difference  between  the  true  and  recovered  model  falls 
within  the  estimated  standard  deviation. 


Figure  7:  The  observed  and  predicted  decay  curves  for  four  stations  in  the  synthetic  data  set  inversion.  The 
predicted  vertical  component  of  the  response,  represented  by  the  solid  lines,  are  a  good  match  to  the  artificially 
generated  noisy  data  set. 
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Figure  9:  The  observed  and  predicted  decay  curves  for  four  stations  in  the  105  mm  projectile  UXO  field  data  set 
inversion.  The  predicted  decay  of  the  vertical  component  of  the  measured  voltages  are  represented  by  the  solid 
lines,  and  the  symbols  represent  the  Geonics  EM63  field  measurements. 
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(a)  Location  and  Orientation 


(b)  Decay  Parameters 


mo 

mrec 

Expected 

Parameters 

a"’ 

Northing  (m) 

1.9 

1.77 

1.83 

0.008 

Easting  (m) 

2.15 

2.04 

2.00 

0.005 

Burial 
Depth  (m) 

0.6 

0.47 

0.44 

0.01 

(j)  (degrees) 

45 

10.1 

~  0 

0.9 

6  (degrees) 

45 

84.7 

-  90 

0.13 

rrii 

mo 

rCirec 

a'" 

ki 

7.07 

76.8 

2.5 

Pi 

1.00 

0.74 

0.04 

7i 

3.16 

31.8 

8.3 

k2 

7.07 

29.2 

3.7 

P2 

1.00 

1.08 

0.07 

72 

3.16 

6.1 

1.3 

(c)  Diagnostics 


Diagnostic 

Result 

Conclusion 

p 

0.91 

permeable 

ki/k2 

2.63 

rod-like 

P1IP2 

0.69 

rod-like 

Table  2:  Recovered  parameters  for  the  field  data  inversion.  Table  (a)  demonstrates  that  the  inversion  was 
successful  in  obtaining  the  approximate  location  and  orientation  of  the  target.  Table  (b)  lists  the  recovered  decay 
parameters  of  the  two  dipoles.  Table  (c)  lists  the  results  of  applying  the  identification  diagnostics  to  the  recovered 
decay  parameters.  Application  of  the  diagnostics  indicates  that  the  buried  target  is  permeable  and  rod”like  and 
therefore  a  candidate  for  UXO. 
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1  EM63  Data  Collection  at  the  ERDC  UXO  Test  Site 


A  series  of  TEM  measurements  was  carried  out  between  April  5  and  April  21,  2000  at  the  Engineering 
Research  Development  Centre  in  Vicksburg,  Mississippi  to  evaluate  the  applicability  and  performance  of 
the  forward  modelling  and  inversion  algorithm  outlined  in  Part  I  of  this  report.  The  Geonics  EM63  TEM 
sensor  was  used  for  all  the  data  collection  in  this  investigation.  The  EM63  is  a  multi-time  channel  time 
domain  unit  consisting  of  a  Im  x  Im  square  transmitter  coil  and  a  co-axial,  horizontal  50cm  diameter 
receiver  coil  mounted  on  a  two-wheel  trailer.  More  information  about  the  EM63  sensor  can  be  found  in 
“EM63  Full  Time  Domain  Electromagnetic  UXO  Detector:  Operating  Instructions”  (2000) . 


Plan  Measurements 

A  series  of  surveys  were  conducted  over  single  targets  seeded  at  the  WES  UXO  test  site  (Figure  11).  A 
4m  X  4m  square  centered  at  (18m  N,  24m  E)  was  chosen  for  the  surveys.  Prior  to  seeding  the  individual 
targets,  EM63  and  EM61-HH  surveys  were  carried  out  to  ensure  that  the  area  was  “quiet”,  i.e.  to  ensure 
the  area  did  not  contain  metallic  scrap.  The  borders  of  the  4m  x  4m  square  were  marked  with  string  to 
indicate  the  extent  of  the  survey  area.  Lines  for  each  survey  were  run  in  a  AT  -  5  direction  with  a  line 
spacing  of  50cm.  The  location  of  the  sensor  was  measured  more  accurately  by  marking  survey  lines  at 
Im  spacing  with  string  and  by  dropping  a  plumb  line  from  the  center  of  the  receiver/transmitter  loop 
pair.  The  EM63  was  set  to  record  a  time  decay  curve  at  10cm  intervals  triggered  by  the  odometer  in  the 
EM63  trailer  wheels. 


Targets  measured  in  this  survey  setup  included  several  UXO  (37  mm  projectile,  60  mm  mortar,  81 
mm  mortar,  stokes  mortar,  105  mm  mortar  projectile,  155  mm  projectile)  and  a  variety  of  scrap  items 
excavated  during  a  UXO  remediation  project  at  Camp  Croft,  Maryland.  These  targets  were  placed  at 
approximately  the  center  of  the  grid  {2mN,2mE),  and  at  depths  up  to  75cm.  In  all  cases  the  strike  of 
the  target  was  parallel  to  the  line  direction.  In  order  to  save  time,  the  soil,  removed  when  digging  a 
hole  for  the  target,  was  not  replaced  over  the  target.  A  wooden  plank  was  placed  over  the  hole  in  which 
the  object  was  laid.  Targets  were  generally  measured  in  three  orientations:  horizontal,  vertical,  and  an 
intermediate  angle. 


Decay  Measurements 

A  controlled  set  of  experiments  was  carried  out  to  carefully  examine  how  the  secondary  field  of  a  target 
decays  as  a  function  of  depth  and  orientation.  For  these  experiments  we  required  measurements  with  the 
EM63  transmitter /receiver  coil  directly  above  the  center  of  various  targets  positioned  at  several  depths 
and  orientations. 


In  order  to  accurately  and  quickly  position  each  target,  Jose  Llopis  at  WES  designed  and  built  a  target 
holder  (Figure  12(a)).  The  jig  was  made  of  wood  and  glue,  contained  no  metallic  materials,  and  could 
orient  each  target  at  15  degree  increments  from  vertical  to  horizontal.  The  size  of  the  jig  made  it  difficult 
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to  bury.  Therefore  the  jig  was  only  partially  buried  and  2x6  planks  were  used  to  adjust  the  height 
between  the  EM63  and  the  target  (Figure  12(b)).  The  height  of  the  planks  was  varied  by  changing  the 
number  of  cinder  blocks  used  to  elevate  the  planks.  The  majority  of  measurements  were  taken  at  two 
sensor  heights:  Z  ~  50cm  and  Z  ~  100cm  from  the  center  of  the  receiver  to  the  center  of  the  target, 
where  Z  is  the  vertical  distance  from  the  center  of  the  target  to  the  center  of  the  receiver  loop.  Additional 
measurements  with  Z  ~  75cm  were  made  on  a  subset  of  the  targets.  TEM  soundings  were  recorded  for 
several  UXO  items  ranging  from  a  37mm  projectile  to  a  105  mm,  as  well  as  several  scrap  items.  Samples 
soundings  are  plotted  in  Figure  13. 


Several  of  the  smaller,  non-UXO  items  did  not  produce  a  significant  response  when  placed  in  the  jig. 
Therefore,  they  were  measured  using  the  setup  of  Figure  14(a).  Each  target  was  placed  in  two  orientations 
and  at  ground  level.  To  facilitate  the  collection  of  cleaner  data,  a  pair  of  6x2  inch  planks  was  placed 
on  the  ground  to  provide  a  level  surface  for  the  EM63  to  be  pulled  along.  The  measured  voltage  curves 
for  scrap  targets  1  to  8  are  plotted  in  Figures  14(b)  and  (c).  The  line  profiles  for  the  first  time  channel 
{t  =  180/is)  and  the  tenth  time  channel  {t  =  0.72msec)  are  plotted  in  Figure  15. 


EM63  TEM  Survey  of  the  ERDC  UXO  Test  Site 

The  ERDC  UXO  Test  Site  consists  of  a  30m  x  95m  grassy  field  containing  seeded  UXO  and  scrap.  On 
the  eastern  end  of  the  site  there  are  approximately  25  targets  buried  along  4  north-south  oriented  rows. 
An  EM63  survey  was  carried  out  between  75mE  and  93.5mE.  The  survey  was  completed  with  Im  line 
spacing  and  10cm  station  spacing.  In  order  to  maintain  straight  lines,  string  was  laid  at  2m  intervals. 
Two  data  sets  were  obtained.  The  first  data  set  was  from  IbmE  to  86.5mE,  and  the  second  data  set 
was  from  82mE  to  93.5mE.  Figures  16(a)  and  16(b)  contain  contour  plots  for  the  first  (0.18ms)  and 
tenth  (0.72ms)  time  channels.  The  plots  were  obtained  by  stitching  together  the  two  data  sets.  This  was 
achieved  by  averaging  the  data  collected  on  duplicate  stations.  For  presentation  purposes,  the  data  in 
each  contour  plot  was  thresholded  from  below  at  a  ImV  level  and  truncated  above  at  levels  indicated  in 
the  plot  titles. 


EM63  TEM  Survey  of  the  ERDC  UXO  Test  Site 
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Figure  11:  4m  x  4m  area  on  which  a  series  of  EM63  surveys  were  taken  over  different  targets  at  several  depths 
and  orientations. 


Figure  12:  (a)  Photograph  of  the  wooden  target  holder  for  EM63  measurements.  A  105mm  projectile  is  standing 
beside  the  jig.  (b)  Measurement  procedure  when  using  the  jig. 
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Figure  13:  (a)  Measurement  of  a  60  mm  and  81  mm  in  jig.  (b)  Measurement  of  Scrap  5  and  a  Steel  Disk.  Unfilled 
symbols  (e.g.  ’o’)  indicate  negative  data  that  has  its  absolute  value  plotted.  Predicted  response  is  obtained  by 
fitting  the  measurement  with  Equation  24. 


EM63  TEM  Survey  of  the  ERDC  UXO  Test  Site 
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Figure  14:  (a)  Photograph  of  setup  used  to  measure  response  of  smaller  scrap  targets,  (b)  Voltage  decay  mea¬ 
surements  for  scrap  targets  1  to  4.  (c)  Voltage  decay  measurements  for  scrap  targets  5  to  8.  Unfilled  symbols 
(e.g.  ’o’)  indicate  negative  data  that  has  its  absolute  value  plotted. 
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station  (cm) 


Station  (cm) 


Figure  15:  Measured  response  of  scrap  targets  1  to  8  along  a  line.  Panel  (a)  and  (b)  plot  the  response  along  a 
line  for  the  first  time  channel  (t  =  0.18msec).  Panel (c)  and  (d)  plot  the  response  along  a  line  for  the  tenth  time 
channel  {t  =  0.72msec).  Each  target  is  located  at  station  0cm. 


EM63  TEM  Survey  of  the  ERDC  UXO  Test  Site 
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Truncate  =  100 


2  Verification  of  the  Dipole  Model 


In  Part  I  of  this  report  an  approximate  forward  model  was  proposed  for  the  TEM  response  of  a  buried 
axi-symmetric  metallic  target.  In  this  forward  model  the  secondary  field  is  approximated  by  a  pair  of 
orthogonal  and  independently  decaying  dipoles,  whose  strengths  are  proportional  to  the  projection  of 
the  primary  field  onto  their  direction.  The  decay  of  each  dipole  moment  is  governed  by  the  magnetic 
polarization  tensor  M.  The  magnetic  polarizability  tensor  is  independent  of  transmitter/receiver/target 
geometry  and  is  a  function  of  the  physical  characteristics  of  the  target  alone.  In  our  previous  work  we 
outlined  a  technique  for  exploiting  M  as  a  tool  for  characterizing  the  shape  of  target  as  either  rod-like  or 
plate-like,  and  if  the  target  is  ferrous  or  non-ferrous  (Pasion  (1999)).  These  results  can  only  be  applied 
with  confidence  for  measurement  configurations  where  the  forward  modelling  is  applicable. 

This  section  contains  a  series  of  tests  designed  to  verify  that  the  magnetic  polarization  tensor  is  indeed 
independent  of  transmitter/receiver/target  geometry.  Each  test  follows  the  same  procedure.  First,  the 
target  is  measured  in  two  orientations:  with  the  axis  of  symmetry  parallel  to  the  primary  field  and  then 
perpendicular  to  the  primary  field.  These  measurements  allow  us  to  extract  the  target’s  two  characteristic 
decay  curves  Li  {t)  and  L2  {t)  that  define  the  magnetic  polarization  tensor.  Equipped  with  M,  we  can 
then  predict  the  TEM  response  for  various  locations  and  orientations  and  compare  these  results  with 
measured  data.  In  this  section  we  evaluate  the  accuracy  of  these  predictions. 


Obtaining  the  Polarizability  Tensor 

Constructing  the  polarizability  tensor  requires  obtaining  the  dipole  decay  functions  Li  (t)  and  L2  (t). 
The  decay  functions  can  be  isolated  by  making  two  measurements:  (1)  with  the  primary  field  parallel 
to  the  z'  axis  of  symmetry,  and  (2)  with  the  primary  field  perpendicular  to  the  z'  axis  of  symmetry. 


Figure  17  illustrates  the  arrangement  of  the  EM63  and  target  that  we  used  to  obtain  the  decay  functions. 
In  this  geometry,  the  measured  voltage  is  then 

V  (t)  =  2  [Li  {t)  cos®  e  +  L2  (t)  sin®  0]  (24) 

where  Z  is  the  distance  between  the  center  of  the  receiver  loop  and  the  center  of  the  target,  and  k  is  a 
constant  that  depends  on  the  size  of  the  receiver  and  transmitter  loops,  the  number  of  turns  in  each  loop, 
and  the  transmitter  current. 


When  the  target’s  z'  axis  of  symmetry  is  parallel  to  the  inducing  field  (0  —  0  degrees.  Figure  17(a)),  only 
the  mi  (i)  dipole  is  excited,  and  the  measured  voltage  is 


V”  (t)  =  (f)  = 


2k 


B^{Z) 


fci(t  +  ai)  e‘ 


*/7i 


(25) 


and  when  the  target’s  z'  axis  of  symmetry  is  perpendicular  to  the  inducing  field  (0  —  90  degrees,  Fig- 
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ure  17(b)),  then  only  the  m2  (t)  dipole  is  excited  and  the  observed  voltage  is 


yJ-  (t)  =  it) 


2k 


B^{Z) 


k2  {t  +  0:2) 


-02  g-^/72 


(26) 


The  voltage  curves  recorded  by  the  EM63  for  the  parallel  and  perpendicular  responses  can  be  obtained 
by  using  a  scaled  down  version  of  the  inversion  algorithm  to  recover  the  decay  parameters  (Pasion,  1999). 
Figure  17(c)  and  (d)  give  the  result  of  this  procedure  applied  to  a  105  mm  projectile  and  scrap  14.  As 
expected  for  a  rod-like  target  the  Li  component  of  the  polarization  tensor  has  a  greater  magnitude  than 
the  L2  component.  Figure  18  has  the  recovered  Li  (t)  and  L2  {t)  for  the  targets  placed  in  the  jig.  Several 
of  the  curves  have  been  extended  by  a  “dash-dot”  line  that  indicates  an  extrapolated  portion  of  Li  (t) 
and  1/2  (f )  • 


The  accuracy  of  this  procedure  to  obtain  Li  (t)  and  L2  {t)  will  depend  on  experimental  error  and  model 
error.  Possible  experimental  errors  include:  (1)  inaccurate  measurement  of  Z]  (2)  inaccurate  placement 
of  the  target  beneath  the  receiver  loop,  i.e.  the  center  of  the  target  must  be  placed  along  the  vertical 
axis  passing  through  the  center  of  the  loop;  (3)  tilting  of  the  EM63  trailer  such  that  the  primary  field 
isn’t  vertical.  Modelling  errors  describe  instances  where  assumptions  of  the  forward  model  are  violated, 
such  as:  (1)  uniformity  of  the  primary  field  in  the  volume  of  the  target;  (2)  representing  the  response  as 
a  point  dipole;  and  (3)  absence  of  fore-aft  symmetry. 


Test  1:  Reproducing  the  time  sounding  at  intermediate  orientation 


Once  we  determine  the  parameters  of  the  magnetic  polarization  tensor,  we  can  forward  model  the  param¬ 
eters  to  obtain  the  TEM  response  for  any  location,  depth,  and  orientation.  In  this  test  we  investigate  how 
accurately  we  could  predict  the  decay  of  the  secondary  field  at  an  arbitrary  orientation  for  a  recovered 
magnetic  polarization  tensor. 


Equation  24  describes  the  measured  voltage  in  a  receiver  coil  directly  above  a  target  illuminated  by  a 
purely  vertical  primary  field.  Using  Equation  24  and  the  decay  functions  Li  {t)  and  L2  (t)  obtained  in  the 
analysis  of  the  previous  section,  we  predict  the  voltage  response  and  compare  it  to  measurements  using 
the  EM63.  The  set  of  measurements  we  use  for  comparison  were  those  obtained  using  the  target  holder. 
Figures  19  and  20  compare  the  measured  responses  at  different  angles  and  the  response  predicted  by  the 
forward  model  for  an  81mm  mortar  without  fins.  Figure  19  has  the  measured  voltage  curves  at  different 
heights  from  the  sensor.  Li  (t)  and  L2  obtained  by  fitting  these  curves  with  equations  25  and  26, 
are  then  forward  modelled  and  plotted  in  Figure  20.  The  top  two  panels  of  Figure  20  demonstrate  the 
procedure  on  an  81mm  mortar  (without  fins)  located  approximately  55cm  beneath  the  receiver  loop.  At 
this  distance  (equivalent  to  approximately  10  cm  below  the  surface)  we  see  that  the  representation  is 
only  moderately  good  at  reproducing  the  data  at  the  different  angles.  When  we  repeat  the  procedure 
for  data  collected  approximately  100cm  beneath  the  receiver,  we  see  that  the  model  does  a  better  job  of 
predicting  the  data.  This  is  not  surprising,  since  we  would  expect  the  modelling  of  a  compact  metallic 
object  as  a  dipole  to  become  more  applicable  as  we  move  further  from  the  source/receiver  loop. 

Test  1:  Reproducing  the  time  sounding  at  intermediate  orientation 
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Test  2:  Reproducing  the  spatial  behavior 

In  our  second  test  we  focus  on  how  accurately  we  could  predict  the  spatial  response  of  a  target  using 
the  dipole  model.  For  this  investigation  we  first  use  a  target’s  magnetic  polarization  tensor,  obtained 
in  the  manner  described  in  the  previous  section,  to  predict  the  secondary  field  over  a  survey  line  that 
passes  directly  over  the  target.  In  each  case  the  survey  line  is  co-aligned  with  the  target.  This  predicted 
response  is  then  compared  with  the  measured  response. 

In  Figures  22  to  30,  the  predicted  and  measured  responses  along  a  line  are  plotted  for  several  targets 
in  three  orientations:  vertical,  horizontal,  and  at  an  intermediate  dip  angle.  At  intermediate  angles  the 
target  dips  towards  the  end  of  the  line  (i.e.  dips  downwards  to  the  right  in  Figures  22  to  30).  Figures  22 
to  30  indicate  the  model  is  successful  in  predicting  the  response  along  each  survey  line. 


Test  2:  Reproducing  the  spatial  behavior 
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Time  (ms)  Time  (ms) 


Figure  17:  (a)  EM63  and  target  geometry  used  to  obtain  Li  (t).  (b)  EM63  and  target  geometry  used  to  obtain 
^2  (t).  (c)  Li{t)  and  L2{t)  curves  for  a  105mm  projectile,  (d)  Li{t)  and  I/2  (t)  curves  for  Scrap  14  (See 
Appendix  A  for  a  photo,  size,  and  weight  of  Scrap  14). 
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0.180  1.0  10.0  25.14 


Time  (ms) 

Figure  19:  Measured  voltage  curves  for  a  81  mm  Mortar  without  fins  at  different  heights  from  the  sensor.  Unfilled 
symbols  (e.g.  ’□’)  indicate  negative  data  that  has  its  absolute  value  plotted.  These  curves  were  inverted  to  obtain 
Li  (t)  and  L2  (t).  Clearly,  the  decay  nature  of  the  signal  changes  when  the  mortar  is  brought  close  to  the  sensor. 
The  resulting  Li  (t)  and  L2  (i)  curves  are  used  to  predict  the  decay  at  intermediate  angles  in  Figure  20. 
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Z  55  cm 


Z  -  55  cm 


Z- 100  cm  Z-- 100  cm 


Time  (ms)  Time  (ms) 

Figure  20:  Predicting  the  measured  voltage  decay  curves  for  an  81  mm  mortar  without  fins  at  different  orientations 
e.  Unfilled  symbols  (e.g.  ’□’)  indicate  negative  data  that  has  its  absolute  value  plotted.  The  Li  (t)  and  L2  (t) 
decay  curves  recovered  from  measurements  with  Z  =  60cm  and  Z  =  52cm  (see  figure  19,  respectively,  are  used  to 
predict  the  voltage  curves  in  the  top  two  panels.  The  Li  (t)  and  L2  (t)  decay  curves  recovered  from  measurements 
^  —  102cm  and  Z  —  95cm,  respectively,  are  used  to  predict  the  voltage  curves  in  the  bottom  two  panels. 
The  prediction  of  the  measured  curves  is  more  successful  when  the  target  is  farther  away  from  the  sensors. 


Test  2:  Reproducing  the  spatial  behavior 
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Z  ~  60  cm  Z  ~  60  cm 


Figure  21:  Predicting  the  measured  voltage  decay  curve  for  Scrap  5  using  the  dipole  model. 
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Figure  22:  60mm  Mortar:  6  =0  degrees  (vertical),  Z  =  73cm. 
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Figure  23:  60mm  Mortar:  9  =90  degrees  (horizontal),  Z  =  73cm. 
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Figure  24:  60mm  Mortar:  6  =40  degrees,  Z  =  75cm. 
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Figure  25:  81mm  Mortar:  0  =0  degrees  (vertical),  Z  =  73cm. 
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Figure  26:  81mm  Mortar:  d  =90  degrees  (horizontal),  Z  =  73cm. 
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Figure  27:  81mm  Mortar:  6  =43  degrees,  Z  =  74cm. 
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Figure  28:  Stokes  Mortar:  B  =90  degrees  (vertical),  Z  =  73cm. 
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Figure  29:  Stokes  Mortar:  6  =90  degrees  (horizontal),  Z  =  73cm 
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Figure  30:  Stokes  Mortsu":  9  =54  degrees,  Z  =  73cm. 
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3  Inversion  of  TEM  D3t3  Collected  3t  the  ERDC  Test 
Site 


In  Part  I  of  this  report  we  described  a  method  of  inverting  TEM  data  for  the  13  parameters  of  the 
Two-Dipole  model.  We  formulated  a  non-linear  least  squares  problem  that  involved  minimizing 


#(m)  =  i||Wd 


(27) 


where  F[m]  is  the  forward  modelled  data,  d^’'’®  is  the  observed  data,  and  $  is  the  least  squares  objective 
function  that  measures  how  closely  our  predicted  data  matches  the  observed  data.  Wd  is  the  data 
weighting  matrix 


(Wd)a  = 


1 

pdi  -h  e 


whore  p  is  typically  a  percentage  and  e  is  a  constant  that  characterizes  ambient  noise. 


(28) 


In  Part  I  we  performed  two  transformations  to  the  model:  (1)  a  simple  linear  scaling  to  the  model  m 
to  ensure  that  each  component  of  m  is  approximately  unity,  and  (2)  replacing  each  of  the  time  decay 
parameters  mi  by  a  squared  variable  Wi  such  that  mi  —  xj  to  ensure  positivity  of  the  decay  parameters. 
For  the  inversions  in  Part  II  of  this  report  we  retain  the  linear  scaling  of  the  model,  but  we  no  longer  solve 
the  square- variable  unconstrained  problem.  Instead  we  reformulate  the  inversion  so  that  it  can  be  solved 
by  a  constrained  optimization  procedure.  We  use  a  Projected  BFGS  algorithm,  described  in  ‘Iterative 
Methods  for  Optimization”  (Kelley,  1999),  to  minimize  the  objective  function.  In  this  algorithm  we  are 
required  to  supply  a  starting  model  as  well  as  upper  and  lower  constraints  for  each  variable. 


Selection  of  Starting  Parameters 

The  first  step  of  this  inversion  procedure  is  to  make  a  starting  guess  for  the  13  model  parameters: 
m  =  [X,y,Z,0,^,A:i,ai,/?i,7i,A;2,a2,^2,72] 

The  success  of  a  local  inversion  procedure,  as  well  as  the  rate  of  convergence  towards  a  solution  is 
dependent  on  the  quality  of  the  initial  guess  (Bard,  1974).  The  following  section  will  discuss  a  number 
of  simple  data  preprocessing  strategies  that  enables  us  to  make  a  reasonable  initial  guess  of  the  above 
parameters. 

Location  on  Survey  (X,  Y) 

A  natural  idea  for  determining  the  location  for  a  target  would  be  to  find  the  location  at  which  the 
maximum  signal  occurs.  This  technique  would  produce  the  exact  location  for  a  spherical  target,  where 
the  plot  of  the  vertical  (z)  secondary  field  produces  a  circular  ”bulls-eye”  anomaly  with  a  maximum 
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located  directly  over  the  target  (Figure  31(a)).  However  rod-like  and  plate-like  targets  can  produce  non- 
symmetric  anomalies  with  multiple  local  maxima  (Figure  31(b)),  In  these  cases,  search  for  the  maximum 
signal  will  give  an  inaccurate  estimation  of  target  location.  To  account  for  these  possible  anomaly  types, 
we  use  an  approximation  of  the  first  spatial  moment  of  the  data  to  estimate  the  location  on 

the  survey: 

^  and  =  EL  Vi  (ti)  yi 

EZiVi(ti)  ZZiVidi) 

where  there  are  N  total  stations  in  the  survey  and  Vi  (h)  is  the  voltage  measured  at  the  i*'*  station  at 
location  (xi,j/i),  and  at  the  first  time  channel.  Figure  31  compares  this  technique  to  simply  using  the 
maximum  value. 


Depth  of  target  (Z) 

The  full  width  at  half  maximum  (FWHM)  is  a  standard  technique  for  estimating  a  target’s  depth.  In 
general,  the  plotted  anomaly  of  non-spherical  targets  will  not  be  circularly  symmetric.  In  these  cases,  it 
is  not  clear  which  profile  to  use  for  calculating  the  FWHM.  Therefore  we  modify  the  FWHM  technique 
by  plotting  points  of  the  secondary  field  anomaly  at  an  early  time  channel  that  are  greater  than  one-half 
of  the  maximum  signal.  The  estimate  for  the  depth  of  the  target  is  then 


Estimated  Depth  = 


(31) 


where  A  is  the  area  of  the  portion  of  the  secondary  field  anomaly  greater  than  one-half  of  the  maximum 
signal. 


Orientation  (</>,  0) 

The  orientation  angles  <j)  and  9  are  not  estimated  prior  to  the  inversion.  The  starting  values  for  (/>  and  9 
set  to  45  degrees. 


Decay  Parameters  {ki,ai,Pi,ji) 

The  decay  parameters  contained  within  the  decay  functions  Li  (t)  and  L2  (t)  are  estimated  in  the  same 
manner  described  in  Pasion  (1999).  The  decay  parameters  are  estimated  by  taking  the  estimate  of  the 
target  location  and  extracting  the  decay  curve  at  that  location.  The  decay  curve  is  then  fit  with  the 
function 

^  (f)  =  (d  -H  tf  e‘/^  (32) 

where  is  the  estimated  depth.  The  starting  decay  parameters  for  both  dipoles  1  and  2  are  then 
(k,a,P,j).  Equation  32  is  obtained  by  setting  Li  =  L2  in  equation  24. 

Selection  of  Upper  and  Lower  Limits 

The  non-linear  least  squares  objective  function  (equation  27)  is  minimized  using  an  optimization  algorithm 
that  incorporates  box  constraints.  Box  constraints  are  the  lower  and  upper  limits  that  we  place  on  the 
model  parameters  that  represent  the  minimum  and  maximum  possible  values  of  the  model  parameters. 
The  upper  and  lower  limits  used  for  the  model  parameters  are  summarized  in  Table  3. 

Selection  of  Upper  and  Lower  Limits 


53 


Location  {X,  Y)  and  Depth  {Z) 

For  the  inversion  results  presented  here  conservative  constraints  were  applied.  The  location  was  con¬ 
strained  to  a  Im  square  box  centered  at  {X,y)  =  {2m,  2m),  and  the  depth  was  constrained  between 
Ziow  =  0.43m  which  is  the  level  of  the  ground,  and  Zhigh  =  l-Sm. 


Orientation  (0, 9) 

Since  we  do  not  estimate  the  orientation  of  the  target  based  on  the  plotted  anomaly,  the  upper  and  lower 
limits  of  <t>  and  e  are  determined  to  be  the  limits  that  allow  for  all  possible  orientations  of  the  target. 
Therefore,  -90  <<t>  <9Q  and  0  <  ^  <  180. 

Decay  Parameters  {ki,ai,  I3i,^i) 

In  verifying  the  applicability  of  the  approximate  forward  model  (Section  2),  the  decay  parameters  were 
recovered  for  different  scrap  and  UXO  targets.  From  this  analysis,  and  the  analysis  reported  in  Pasion 
(1999),  upper  and  lower  limits  for  the  decay  parameters  were  chosen.  The  limits  are  listed  in  Table  3(b). 

Application  to  Field  Data  Sets 

In  this  section  we  present  the  results  of  applying  the  inversion  algorithm  to  five  targets:  155  mm,  105 
mm,  81  mm,  60  mm,  and  Stokes  mortar.  A  list  of  inversion  results  presented  in  this  report  can  be  found 
in  Table  4,  where  the  target  type,  depth,  and  orientation  for  each  data  set  inverted  in  this  section  are 
given.  Each  inverted  data  set  consisted  of  soundings  collected  on  a  2  m  x  2  m  grid  centered  on  the  target, 
containing  5  lines  running  North-South  separated  at  50  cm  line  spacing,  with  stations  located  at  5  cm 
intervals  along  each  line.  Inversion  results  can  be  found  in  Figures  32  to  48.  In  Panel  (a)  of  each  figure  is 
a  plan  view  comparison  of  the  observed  and  predicted  data  for  time  channels  1,  8,  15,  and  23.  Panel  (b) 
compares  the  predicted  and  observed  decay  at  three  stations:  Line  1  m.  Station  Im;  Line  2  m.  Station  2 
m;  and  Line  2.5  m.  Station  2.5  m.  The  recovered  location  and  orientation  are  listed  in  Panel  (c),  and  the 
recovered  decay  parameters  and  diagnostics  applied  to  those  parameters  are  reported  in  Panel  (d).  The 
algorithm  was  successful  in  identifying  each  target  as  rod-like  and  permeable  for  all  the  cases  presented 
here. 
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Easting  (m)  Easting  (m) 


Figure  31:  Finding  the  location  of  a  (a)  sphere  and  (b)  horizontal  stokes  mortar.  The  location  obtained  by 
using  the  maximum  of  the  signal  is  indicated  by  the  star,  and  the  location  obtained  by  using  the  “center-of 
mass”  is  indicated  by  the  triangle.  The  center-of-mass  is  much  better  at  determining  the  location  of  the  target 
when  the  plotted  anomaly  does  not  have  the  appearance  of  a  “bulls-eye”  (e.g.  a  sphere).  Panels  (c)  and  (d) 
plot  those  data  points  that  are  greater  than  half  the  maximum  signal.  Panel  (d)  indicates  that  the  traditional 
“Full-Width-at-Half-Maximum”  would  be  difficult  to  utilize  for  non- “bulls-eye”  anomalies. 


(a)  Location  and  Orientation 


(b)  Decay  Parameters 


rrii 

ki 

0.1 

200 

O^i 

0.001 

0.02 

Pi 

0.3 

1.5 

li 

2 

30 

rrii 

rrii 

X  (m) 

1.5 

2.5 

Y(m) 

1.5 

2.5 

Z  (m) 

0.43 

1.5 

<f)  (degrees) 

-90 

90 

9  (degrees) 

0 

180 

Table  3:  Upper  and  lower  limits  used  in  inversion  examples  in  Part  II  of  this  report. 
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Target 

Z  (cm) 

9  (degrees) 

Figure  No. 

155  mm 

98 

0 

32 

155  mm 

97 

52 

33 

155  mm 

91 

90 

34 

105  mm 

86 

0 

35 

105  mm 

87 

90 

36 

105  mm 

89 

59 

37 

105  mm 

89 

40 

38 

Stokes  Mortar 

73 

54 

39 

Stokes  Mortar 

73 

0 

40 

Stokes  Mortar 

73 

90 

41 

81  mm 

73  1 

0 

42 

81  mm 

74 

43 

43 

81  mm 

73 

90 

44 

60  mm 

58 

57 

45 

60  mm 

58 

90 

46 

60  mm 

59 

0  (nose  up) 

47 

60  mm 

59 

0  (nose  down) 

48 

Table  4:  Targets 
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Figure  32:  Recovered  parameters  for  tlie  field  data  inversion  of  a  155mm  projectile. 


Location  and  Orientation  (d)  Decay  Parameters  and  Diagnostics 
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Figure  33:  Recovered  parameters  for  the  field  data  inversion  of  a  155mm  projectile. 


Figure  36:  Recovered  parameters  for  the  field  data  inversion  of  a  105mm  projectile. 


Observed  Predicted 


Figure  37:  Recovered  parameters  for  the  field  data  inversion  of  a  105mm  projectile. 


Figure  39:  Recovered  parameters  for  the  field  data  inversion  of  a  stokes  mortar. 


Stokes  Mortar:  0  =  0  degrees  (vertical),  Z  =  73cm 
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Location  and  Orientation 


Application  to  Field  Data  Sets 


Figure  43:  Recovered  parameters  for  the  field  data  inversion  of  an  81  mm  mortar  without  fins, 


81  mm;  0  =  90  degrees  (horizontal),  Z  =  73cm 
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44:  Recovered  parameters  for  the  field  data  inversion  of  an  81  mm  mortar  without  fins. 


57  deg 


Application  to  Field  Data  Sets 


Figure  45:  Recovered  parameters  for  the  field  data  inversion  of  a  60mm  mortar. 


60  mm:  9  =  90  degrees  (horizontal),  Z  =  58cm 


Application  to  Field  Data  Sets 


Figure  46:  Recovered  parameters  for  the  field  data  inversion  of  a  60mm  mortar. 


60  mm:  6  =  0  degrees  (vertical,  nose  up),  Z  =  bOcm 


Application  to  Field  Data  Sets 


Figure  47:  Recovered  paxameters  for  the  field  data  inversion  of  a  60mm  mortar 
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4  References 


A  Target  Descriptions 


Target 

Mass 

(kg) 

Length 

(cm) 

Width 

/Thickness  (cm) 

155  mm 

60 

15.5 

105  mm 

14.05 

21 

10.5 

Stokes  Mortar 

3.55 

36 

16 

81  mm  (no  fins) 

3.3 

26 

8.1 

Rusted  Mortar 

2.85 

36 

29 

60  mm 

1.352 

26 

6 

37  mm 

0.839 

11.4 

3.6 

Disk 

0.88 

8  (diameter) 

2.2  thick 

Scrap  1 

0.08 

13  (diameter) 

^2  thick 

Scrap  2 

0.026 

6.4 

6.3 

Scrap  3 

0.055 

11 

6.4 

Scrap  4 

0.091 

19.4 

2.5 

Scrap  5 

0.172 

16.5 

3 

Scrap  6 

0.256 

21 

4.8 

Scrap  7 

0.069 

12.5 

2 

Scrap  8 

0.075 

7 

3 

Scrap  9 

19 

9 

Scrap  10 

0.032 

5 

3.5 

Scrap  12 

0.186 

9.6 

4.4 

Rocket  Fins 
(Scrap  13) 

0.83 

21.4 

6.5 

Blown  Mortar 
(Scrap  14) 

0.939 

16 

7 

Scrap  15 

0.431 

16 

5.6 

A  Target  Descriptions 


B  Recovered  Polarization  Dyadic  for  Various  Targets 


The  magnetic  polarization  dyadic  for  the  TEM  problem  is  defined  as 

'L2  (t)  0  0 

M  =  0  Liit)  0 

0  0  Li{t) 

where 

Li  (t)  ^ki{t  + 

L2it)^k2it  +  a2r^^e-*/^^ 

The  recovered  magnetic  polarization  dyadic  acquired  for  a  number  of  targets  are  reported  here.  The 
procedure  followed  to  obtained  these  results  is  described  in  Section  2  of  this  report.  A  photo  and 
description  of  each  target  can  be  found  in  Appendix  A. 


Tjme  (ms) 


Figure  51:  105  mm  Projectile 
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60  mm 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

3.38 

0.79 

a 

0.019 

0.02 

p 

0.90 

1.19 

7 

2.55 

3.59 

Figure  55:  60  mm  Mortar 


(b)  Decay  Parameters 


mi 

Dipole  1 

Dipole  2 

k 

2.33 

0.14 

a 

0.02 

0.001 

13 

0.40 

1.46 

7 

3.49 

2.36 

Figure  56:  37  mm  Projectile 


(b)  Decay  Parameters 
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Figure  57:  Fins  from  a  2.75  inch  Rocket 
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Scrap  14 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

1.78 

0.89 

a 

0.02 

0.02 

0 

0.80 

1.11 

7 

8.99 

2.14 

Figure  58:  Scrap  14 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

1.26 

na 

a 

0.02 

na 

0 

0.81 

na 

7 

2.08 

na 

Figure  59:  Scrap  5.  The  response  of  the  m2  dipole  did  not  generate  a  response  above  the  noise  level  of  the  EM63. 


(b)  Decay  Parameters 


rrii 

Dipole  1 

Dipole  2 

k 

0.23 

1.45 

a 

0.001 

0.02 

0 

1.44 

0.99 

7 

2.57 

2.00 

Figure  60:  Steel  Disk 
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